##########################################
### LANE AND SCHOENHERR CITE AND SWAY? ###
##########################################

# Purpose:
# 	- Replicate Tables B1-B6, B8 in the supplemental appendix

# Data Sources:
#	- Analysis2DataAppendixPart1.csv
#	- Analysis2DataAppendixPart2.csv
#	- Analysis2DataAppendixPart3.csv

# Run on:
#	- R 4.4.0 ("Puppy Cup") on macOS Sequoia 15.1.1.

# Associated Codebook:
#	- Analysis2Codebook.pdf

###########################################
###########################################
###########################################

library(foreign)
library(arm)
library(readxl)
library(lattice)
library(tidyverse)
library(stargazer)
library(faraway)
library(lme4)
library(haven)

#####################################
### STEP 0: SET WORKING DIRECTORY ###
#####################################

setwd()

########################################
### STEP 1: PULL IN THE MODEL 2 DATA ###
########################################

model2data <- read.csv("Analysis2DataAppendixPart1.csv")

############################################
### STEP 2: PRE-PROCESS SOME OF THE DATA ###
############################################

model2data$logPetBriefExperience <- log(model2data$petBriefExperience + 1)
model2data$logRespBriefExperience <- log(model2data$respBriefExperience + 1)

model2data$logPetNumCites <- log(model2data$petNumCites + 1)
model2data$logRespNumCites <- log(model2data$respNumCites + 1)

# check all the variables to make sure no NAs before running:
table(is.na(model2data$voteWithPet))
table(is.na(model2data$petUnnamedCiteCount))
table(is.na(model2data$petMajCallOutCount))
table(is.na(model2data$petNotMajOpinCount))
table(is.na(model2data$petDissent))
table(is.na(model2data$petConcurrence))
table(is.na(model2data$respUnnamedCiteCount))
table(is.na(model2data$respMajCallOutCount))
table(is.na(model2data$respNotMajOpinCount))
table(is.na(model2data$respDissent))
table(is.na(model2data$respConcurrence))
table(is.na(model2data$ideoAlign))
table(is.na(model2data$logPetBriefExperience))
table(is.na(model2data$logRespBriefExperience))
table(is.na(model2data$sgParty))
table(is.na(model2data$helpOSG))
table(is.na(model2data$amiciNet))
table(is.na(model2data$petFormerClerk))
table(is.na(model2data$respFormerClerk))
table(is.na(model2data$netStatus))
table(is.na(model2data$lcDisagreement))
model2data <- model2data %>% filter(!is.na(lcDisagreement))
table(is.na(model2data$CSI))
table(is.na(model2data$pastExpertise))
table(is.na(model2data$petNumCites))
table(is.na(model2data$respNumCites))
table(is.na(model2data$petQuestionsOA))
table(is.na(model2data$respQuestionsOA))
table(is.na(model2data$justice))
table(is.na(model2data$issueArea))
table(is.na(model2data$term))

# convert the three levels to factors
model2data$justice <- as.factor(model2data$justice)
model2data$issueArea <- as.factor(model2data$issueArea)
model2data$term <- as.factor(model2data$term)

#####################
### STEP 3: MODEL ###
#####################

model2data$petUnnamedCiteLog <- log(model2data$petUnnamedCiteCount + 1)
model2data$petMajCallOutLog <- log(model2data$petMajCallOutCount + 1)
model2data$petNotMajOpinLog <- log(model2data$petNotMajOpinCount + 1)
model2data$respUnnamedCiteLog <- log(model2data$respUnnamedCiteCount + 1)
model2data$respMajCallOutLog <- log(model2data$respMajCallOutCount + 1)
model2data$respNotMajOpinLog <- log(model2data$respNotMajOpinCount + 1)
model2data$pastExpertiseLog <- log(model2data$pastExpertise + 1)
model2data$oaQuestDiff <- model2data$petQuestionsOA - model2data$respQuestionsOA
model2data$petExperienceAdvantage <- model2data$petBriefExperience - model2data$respBriefExperience

### TABLE 4 IN MANUSCRIPT ###

model2logCites <- glmer(voteWithPet ~ petUnnamedCiteLog
							+ petMajCallOutLog
							+ petNotMajOpinLog
							+ respUnnamedCiteLog
							+ respMajCallOutLog
							+ respNotMajOpinLog
							+ ideoAlign
							+ pastExpertise
							+ logPetNumCites
							+ logRespNumCites
							+ petExperienceAdvantage
							+ sgParty
							+ lcDisagreement
							+ amiciNet
							+ helpOSG
							+ netStatus
							+ oaQuestDiff
							+ petUnnamedCiteLog * ideoAlign
							+ petMajCallOutLog * ideoAlign
							+ petNotMajOpinLog * ideoAlign
							+ respUnnamedCiteLog * ideoAlign
							+ respMajCallOutLog * ideoAlign
							+ respNotMajOpinLog * ideoAlign
							+ (1 | issueArea)
							+ (1 | term)
							+ (1 | justice),
							data = model2data,
							family = binomial(link = logit),
							glmerControl(optimizer = "Nelder_Mead"))

summary(model2logCites)
display(model2logCites)

####################################
#### STEP 3: MODEL 2 WITH COUNTS ###
####################################

# TABLE B1, COLUMN 2 ###

model2count <- glmer(voteWithPet ~ petUnnamedCiteCount
							+ petMajCallOutCount
							+ petNotMajOpinCount
							+ respUnnamedCiteCount
							+ respMajCallOutCount
							+ respNotMajOpinCount
							+ ideoAlign
							+ pastExpertise
							+ logPetNumCites
							+ logRespNumCites
							+ petExperienceAdvantage
							+ sgParty
							+ lcDisagreement
							+ amiciNet
							+ helpOSG
							+ netStatus
							+ oaQuestDiff						
							+ petUnnamedCiteCount * ideoAlign
							+ petMajCallOutCount * ideoAlign
							+ petNotMajOpinCount * ideoAlign
							+ respUnnamedCiteCount * ideoAlign
							+ respMajCallOutCount * ideoAlign
							+ respNotMajOpinCount * ideoAlign
							+ (1 | issueArea)
							+ (1 | term)
							+ (1 | justice),
							data = model2data,
							family = binomial(link = logit),
							glmerControl(optimizer = "Nelder_Mead"))

summary(model2count)
display(model2count)

####################################
### STEP 3B: WRITE OUT THE TABLE ###
####################################

stargazer(model2logCites, model2count, align = TRUE, omit.stat=c("LL", "ser", "f"))

###################################
### STEP 4: ALT CITE CATEGORIES ###
###################################

# Table B2 in the Supplemental Appendix

#########################################
## STEP 4A: COLLAPSE ACTIVE CITATIONS ###
#########################################

model2data$petAllCallOutCount <- model2data$petMajCallOutCount + model2data$petNotMajOpinCount
model2data$respAllCallOutCount <- model2data$respMajCallOutCount + model2data$respNotMajOpinCount

model2data$petAllCallOutLog <- log(model2data$petAllCallOutCount + 1)
model2data$respAllCallOutLog <- log(model2data$respAllCallOutCount + 1)

### TABLE B2, COLUMN 2 ###

model2allCall <- glmer(voteWithPet ~ petUnnamedCiteLog
							+ petAllCallOutLog
							+ respUnnamedCiteLog
							+ respAllCallOutLog
							+ ideoAlign
							+ pastExpertise
							+ logPetNumCites
							+ logRespNumCites
							+ petExperienceAdvantage
							+ sgParty
							+ lcDisagreement
							+ amiciNet
							+ helpOSG
							+ netStatus
							+ oaQuestDiff
							+ petUnnamedCiteLog * ideoAlign
							+ petAllCallOutLog * ideoAlign
							+ respUnnamedCiteLog * ideoAlign
							+ respAllCallOutLog * ideoAlign
							+ (1 | issueArea)
							+ (1 | term)
							+ (1 | justice),
							data = model2data,
							family = binomial(link = logit),
							glmerControl(optimizer = "Nelder_Mead"))

summary(model2allCall)
display(model2allCall)

#######################################
### STEP 4B: COLLAPSE ALL CITATIONS ###
#######################################

model2data$petAllCitesCount <- model2data$petUnnamedCiteCount + model2data$petMajCallOutCount + model2data$petNotMajOpinCount
model2data$respAllCitesCount <- model2data$respUnnamedCiteCount + model2data$respMajCallOutCount + model2data$respNotMajOpinCount

model2data$petAllCitesLog <- log(model2data$petAllCitesCount + 1)
model2data$respAllCitesLog <- log(model2data$respAllCitesCount + 1)

### TABLE B2, COLUMN 3 ###

model2allCites <- glmer(voteWithPet ~ petAllCitesLog
							+ respAllCitesLog
							+ ideoAlign
							+ pastExpertise
							+ logPetNumCites
							+ logRespNumCites
							+ petExperienceAdvantage
							+ sgParty
							+ lcDisagreement
							+ amiciNet
							+ helpOSG
							+ netStatus
							+ oaQuestDiff
							+ petAllCitesLog * ideoAlign
							+ respAllCitesLog * ideoAlign
							+ (1 | issueArea)
							+ (1 | term)
							+ (1 | justice),
							data = model2data,
							family = binomial(link = logit),
							glmerControl(optimizer = "Nelder_Mead"))

summary(model2allCites)
display(model2allCites)

# note it's a better model to not combine them
# best model is to pull together all the call outs, but that's not a theoretically correct model
# second best model is logged with all three -- that'll be the model we use

###################################################
### STEP 4C: SEPARATE DISSENTS AND CONCURRENCES ###
###################################################

model2data$petDissentLog <- log(model2data$petDissent + 1)
model2data$petConcurrenceLog <- log(model2data$petConcurrence + 1)
model2data$respDissentLog <- log(model2data$respDissent + 1)
model2data$respConcurrenceLog <- log(model2data$respConcurrence + 1)

### TABLE B2, COLUMN 4 ###

model2logCitesSepSplit <- glmer(voteWithPet ~ petUnnamedCiteLog
							+ petMajCallOutLog
							+ petDissentLog
							+ petConcurrenceLog
							+ respUnnamedCiteLog
							+ respMajCallOutLog
							+ respDissentLog
							+ respConcurrenceLog
							+ ideoAlign
							+ pastExpertise
							+ logPetNumCites
							+ logRespNumCites
							+ petExperienceAdvantage
							+ sgParty
							+ lcDisagreement
							+ amiciNet
							+ helpOSG
							+ netStatus
							+ oaQuestDiff
							+ petUnnamedCiteLog * ideoAlign
							+ petMajCallOutLog * ideoAlign
							+ petDissentLog * ideoAlign
							+ petConcurrenceLog * ideoAlign
							+ respUnnamedCiteLog * ideoAlign
							+ respMajCallOutLog * ideoAlign
							+ respDissentLog * ideoAlign
							+ respConcurrenceLog * ideoAlign
							+ (1 | issueArea)
							+ (1 | term)
							+ (1 | justice),
							data = model2data,
							family = binomial(link = logit),
							glmerControl(optimizer = "Nelder_Mead"))

summary(model2logCitesSepSplit)
display(model2logCitesSepSplit)

####################################
### STEP 4D: WRITE OUT THE TABLE ###
####################################

stargazer(model2logCites, model2allCall, model2allCites, model2logCitesSepSplit, align = TRUE, omit.stat=c("LL", "ser", "f"))

####################################
### STEP 5: ISSUE AREA MQ SCORES ###
####################################

# Table B3 in the Supplemental Appendix

# notes:
#	- For Issue Areas 11, 13, 14, just issued the full Court MQ scores due to sparsity
# 	- all other issue area MQ scores calculated using MQ code + 2024 update to MCMCpack
#	- note that issue areas 6, 7, and 12 are not particularly stable
#		- there's less data missing than 11, 13, and 14, but still a lot of issues with stability
#		- mostly driven by the docket drop off in the 1990s

# if there's an NA for model2data$issueMean, drop it out
# means there wasn't enough data to get an estimate
model2dataIssue <- model2data %>% filter(!is.na(issueMean))

##############################
### STEP 5A: RUN THE MODEL ###
##############################

### TABLE B3, COLUMN 2 ###

model2issueMQ <- glmer(voteWithPet ~ petUnnamedCiteLog
							+ petMajCallOutLog
							+ petNotMajOpinLog
							+ respUnnamedCiteLog
							+ respMajCallOutLog
							+ respNotMajOpinLog
							+ ideoAlignIssue
							+ pastExpertise
							+ logPetNumCites
							+ logRespNumCites
							+ petExperienceAdvantage
							+ sgParty
							+ lcDisagreement
							+ amiciNet
							+ helpOSG
							+ netStatus
							+ oaQuestDiff
							+ petUnnamedCiteLog * ideoAlignIssue
							+ petMajCallOutLog * ideoAlignIssue
							+ petNotMajOpinLog * ideoAlignIssue
							+ respUnnamedCiteLog * ideoAlignIssue
							+ respMajCallOutLog * ideoAlignIssue
							+ respNotMajOpinLog * ideoAlignIssue
							+ (1 | issueArea)
							+ (1 | term)
							+ (1 | justice),
							data = model2dataIssue,
							family = binomial(link = logit),
							glmerControl(optimizer = "Nelder_Mead"))

summary(model2issueMQ)
display(model2issueMQ)

####################################
### STEP 5B: WRITE OUT THE TABLE ###
####################################

stargazer(model2logCites, model2issueMQ, align = TRUE, omit.stat=c("LL", "ser", "f"))

##################################################
### STEP 6: KITCHEN SINK MODEL: CLERKS AND CSI ###
##################################################

# Table B4 in the Supplemental Appendix

##############################
### STEP 6A: RUN THE MODEL ###
##############################

### TABLE B4, COLUMN 2 ###
model2big <- glmer(voteWithPet ~ petUnnamedCiteLog
							+ petMajCallOutLog
							+ petNotMajOpinLog
							+ respUnnamedCiteLog
							+ respMajCallOutLog
							+ respNotMajOpinLog
							+ ideoAlign
							+ pastExpertise
							+ logPetNumCites
							+ logRespNumCites
							+ petExperienceAdvantage
							+ sgParty
							+ lcDisagreement
							+ amiciNet
							+ helpOSG
							+ netStatus
							+ oaQuestDiff
							+ petFormerClerk
							+ respFormerClerk
							+ CSI
							+ petUnnamedCiteLog * ideoAlign
							+ petMajCallOutLog * ideoAlign
							+ petNotMajOpinLog * ideoAlign
							+ respUnnamedCiteLog * ideoAlign
							+ respMajCallOutLog * ideoAlign
							+ respNotMajOpinLog * ideoAlign
							+ (1 | issueArea)
							+ (1 | term)
							+ (1 | justice),
							data = model2data,
							family = binomial(link = logit),
							glmerControl(optimizer = "Nelder_Mead"))

summary(model2big)
display(model2big)

####################################
### STEP 6B: WRITE OUT THE TABLE ###
####################################

stargazer(model2logCites, model2big, align = TRUE, omit.stat=c("LL", "ser", "f"))

#################################################
### STEP 7: COSINE SIMILARITY AND READABILITY ###
#################################################

# Table B5 in the Supplemental Appendix

# read in data with Hazelton and Hinkle variables
model2dataRead <- read.csv("Analysis2DataAppendixPart2.csv")

###########################
### STEP 7A: READBILITY ###
###########################

# ONE: Run the full model in the paper, but it's truncated to 2015 because of the data

### TABLE B5, COLUMN 2 ###

model2full2015 <- glmer(voteWithPet ~ petUnnamedCiteLog
							+ petMajCallOutLog
							+ petNotMajOpinLog
							+ respUnnamedCiteLog
							+ respMajCallOutLog
							+ respNotMajOpinLog
							+ ideoAlign
							+ pastExpertise
							+ logPetNumCites
							+ logRespNumCites
							+ petExperienceAdvantage
							+ sgParty
							+ lcDisagreement
							+ amiciNet
							+ helpOSG
							+ netStatus
							+ oaQuestDiff
							+ petUnnamedCiteLog * ideoAlign
							+ petMajCallOutLog * ideoAlign
							+ petNotMajOpinLog * ideoAlign
							+ respUnnamedCiteLog * ideoAlign
							+ respMajCallOutLog * ideoAlign
							+ respNotMajOpinLog * ideoAlign
							+ (1 | issueArea)
							+ (1 | term)
							+ (1 | justice),
							data = model2dataRead,
							family = binomial(link = logit))
							

summary(model2full2015)
display(model2full2015)


### TABLE B5, COLUMN 3 ###

model2readability2015 <- glmer(voteWithPet ~ petUnnamedCiteLog
							+ petMajCallOutLog
							+ petNotMajOpinLog
							+ respUnnamedCiteLog
							+ respMajCallOutLog
							+ respNotMajOpinLog
							+ ideoAlign
							+ pastExpertise
							+ logPetNumCites
							+ logRespNumCites
							+ petExperienceAdvantage
							+ sgParty
							+ lcDisagreement
							+ amiciNet
							+ helpOSG
							+ netStatus
							+ oaQuestDiff							
							+ petReadability
							+ respReadability
							+ petRespCosineSim
							+ petUnnamedCiteLog * ideoAlign
							+ petMajCallOutLog * ideoAlign
							+ petNotMajOpinLog * ideoAlign
							+ respUnnamedCiteLog * ideoAlign
							+ respMajCallOutLog * ideoAlign
							+ respNotMajOpinLog * ideoAlign
							+ (1 | issueArea)
							+ (1 | term)
							+ (1 | justice),
							data = model2dataRead,
							family = binomial(link = logit)
							)

summary(model2readability2015)
display(model2readability2015)

##################################
### STEP 7B: COSINE SIMILARITY ###
##################################

# read in Hazelton and Hinkle cosine similarity data
# this is cosine similarity between party + ACBs
# it's going to drop a ton of stuff out -- not every case has amicus filers on both sides, let alone data for that
model2dataCos <- read.csv("Analysis2DataAppendixPart3.csv")

### TABLE B5, COLUMN 4 ###

model2amicusCosine2015 <- glmer(voteWithPet ~ petUnnamedCiteLog
							+ petMajCallOutLog
							+ petNotMajOpinLog
							+ respUnnamedCiteLog
							+ respMajCallOutLog
							+ respNotMajOpinLog
							+ ideoAlign
							+ pastExpertise
							+ logPetNumCites
							+ logRespNumCites
							+ petExperienceAdvantage
							+ sgParty
							+ lcDisagreement
							+ amiciNet
							+ helpOSG
							+ netStatus
							+ oaQuestDiff
							+ petReadability
							+ respReadability
							+ petRespCosineSim
							+ petAmicusCosine
							+ respAmicusCosine
							+ petUnnamedCiteLog * ideoAlign
							+ petMajCallOutLog * ideoAlign
							+ petNotMajOpinLog * ideoAlign
							+ respUnnamedCiteLog * ideoAlign
							+ respMajCallOutLog * ideoAlign
							+ respNotMajOpinLog * ideoAlign
							+ (1 | issueArea)
							+ (1 | term)
							+ (1 | justice),
							data = model2dataCos,
							family = binomial(link = logit)
							)

summary(model2amicusCosine2015)
display(model2amicusCosine2015)

####################################
### STEP 7C: WRITE OUT THE TABLE ###
####################################

stargazer(model2logCites, model2full2015, model2readability2015, model2amicusCosine2015, align = TRUE, omit.stat=c("LL", "ser", "f"))

############################################
### STEP 8: VARIABLE BREAKDOWN FOR TABLE ###
############################################

# Table B6 in Supplemental Appendix

summary(model2data$voteWithPet)
summary(model2data$petUnnamedCiteLog)
summary(model2data$petMajCallOutLog)
summary(model2data$petNotMajOpinLog)
summary(model2data$respUnnamedCiteLog)
summary(model2data$respMajCallOutLog)
summary(model2data$respNotMajOpinLog)

summary(model2data$ideoAlign)
summary(model2data$pastExpertise)
summary(model2data$logPetNumCites)
summary(model2data$logRespNumCites)

summary(model2data$petExperienceAdvantage)
summary(model2data$sgParty)
summary(model2data$lcDisagreement)

summary(model2data$amiciNet)
summary(model2data$helpOSG)
summary(model2data$netStatus)
summary(model2data$oaQuestDiff)

########################
### STEP 9: MATCHING ###
########################

# Table B8 in the supplemental appendix

# create simplified variables
model2data$petUnnamedPeriod <- ifelse(model2data$petUnnamedCiteCount >= 1, 1, 0)
model2data$petMajCallOutPeriod <- ifelse(model2data$petMajCallOutCount >= 1, 1, 0)
model2data$petNotMajOpinPeriod <- ifelse(model2data$petNotMajOpinCount >= 1, 1, 0)
model2data$petCitePeriod <- ifelse(model2data$petUnnamedPeriod == 1, 1, 
	ifelse(model2data$petMajCallOutPeriod == 1, 1, 
	ifelse(model2data$petNotMajOpinPeriod == 1, 1, 0)))

model2data$respUnnamedPeriod <- ifelse(model2data$respUnnamedCiteCount >= 1, 1, 0)
model2data$respMajCallOutPeriod <- ifelse(model2data$respMajCallOutCount >= 1, 1, 0)
model2data$respNotMajOpinPeriod <- ifelse(model2data$respNotMajOpinCount >= 1, 1, 0)	
model2data$respCitePeriod <- ifelse(model2data$respUnnamedPeriod == 1, 1,
	ifelse(model2data$respMajCallOutPeriod == 1, 1,
	ifelse(model2data$respNotMajOpinPeriod == 1, 1, 0)))

#####################################
### STEP 10A: MATCH ON PETITIONER ###
#####################################

# install the CEM packages
# Code borrowed from Black and Owens (2021) - clerks (per reviewer suggestion)
library(cem)
library(MASS)
library(car)
library(lmtest)

# citing is the treatment variable turn it into a factor
model2data$petCitePeriod <- as.factor(model2data$petCitePeriod)

# Put together the set for CEM
cemData <- model2data %>% select(respCitePeriod, ideoAlign, pastExpertise, petExperienceAdvantage, sgParty, lcDisagreement, amiciNet, helpOSG, netStatus, voteWithPet, petCitePeriod)

# figure out what to drop out for the imbalance calculation
todrop <- c("voteWithPet", "petCitePeriod")

# perform CEM
cemWeights <- cem(treatment = "petCitePeriod", data = cemData, drop = todrop)

cemData2 <- model2data %>% select(respCitePeriod, ideoAlign, pastExpertise, logPetNumCites, logRespNumCites, petExperienceAdvantage, sgParty, lcDisagreement, amiciNet, helpOSG, netStatus, oaQuestDiff, voteWithPet, petCitePeriod, term, issueArea, justice)

# calculate the imbalance pre and post
set.seed(19870302)
pre <- L1.profile(group = cemData$petCitePeriod, data = cemData, drop = todrop, plot = F, M = 500)
post <- L1.profile(group = cemData$petCitePeriod, data = cemData, drop = todrop, plot = F, useCP = pre$CP, weights = cemWeights$w)

pre$medianL1
# 0.5343468

post$medianL1
# 0.12344304

(pre$medianL1-post$medianL1)/pre$medianL1
# 0.769007

################################
### STEP 10B: RUN THE MODELS ###
################################

# do some cleaning
cemData2$issueArea <- as.factor(cemData2$issueArea)
cemData2$justice <- as.factor(cemData2$justice)
cemData2$term <- as.factor(cemData2$term)

# run the models

### TABLE B8, COLUMN 1 ###

model2NoCEM <- glmer(voteWithPet ~ petCitePeriod
							+ respCitePeriod
							+ ideoAlign
							+ pastExpertise
							+ logPetNumCites
							+ logRespNumCites
							+ petExperienceAdvantage
							+ sgParty
							+ lcDisagreement
							+ amiciNet
							+ helpOSG
							+ netStatus
							+ oaQuestDiff
							+ petCitePeriod * ideoAlign
							+ respCitePeriod * ideoAlign
							+ (1 | issueArea)
							+ (1 | term)
							+ (1 | justice),
							data = cemData2,
							family = binomial(link = logit),
							glmerControl(optimizer = "Nelder_Mead"))

summary(model2NoCEM)
display(model2NoCEM)

### TABLE B8, COLUMN 2 ###

model2CEM <- glmer(voteWithPet ~ petCitePeriod
							+ respCitePeriod
							+ ideoAlign
							+ pastExpertise
							+ logPetNumCites
							+ logRespNumCites
							+ petExperienceAdvantage
							+ sgParty
							+ lcDisagreement
							+ amiciNet
							+ helpOSG
							+ netStatus
							+ oaQuestDiff
							+ petCitePeriod * ideoAlign
							+ respCitePeriod * ideoAlign
							+ (1 | issueArea)
							+ (1 | term)
							+ (1 | justice),
							data = cemData2,
							family = binomial(link = logit),
							weights = cemWeights$w,
							glmerControl(optimizer = "Nelder_Mead"))

summary(model2CEM)
display(model2CEM)

######################################
### STEP 10C: WRITE OUT THE TABLES ###
######################################

### TABLE B8 ###
stargazer(model2NoCEM, model2CEM, align = TRUE, omit.stat=c("LL", "ser", "f"))




